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Abstract 

We study the critical behavior of a quenched random-exchange Ising model with competing in- 
teractions on a bcc lattice. This model was introduced in the study of the magnetic behavior of 
Fei-zRuz alloys for ruthenium concentrations x = 0%, x = 4%, x = 6%, and x = 8%. Our study 
is carried out within a Monte Carlo approach, with the aid of a re-weighting multiple histogram 
technique. By means of a finite-size scaling analysis of several thermodynamic quantities, taking 
into account up to the leading irrelevant scaling field term, we find estimates of the critical ex- 
ponents a, /?, 7, and v, and of the critical temperatures of the model. Our results for x = 0% 
are in excellent agreement with those for the three-dimensional pure Ising model in the literature. 
We also show that our critical exponent estimates for the disordered cases are consistent with 
those reported for the transition line between paramagnetic and ferromagnetic phases of both ran- 
domly dilute and ±J Ising models. We compare the behavior of the magnetization as a function 
of temperature with that obtained by Paduani and Branco (2008), qualitatively confirming the 
mean-field result. However, the comparison of the critical temperatures obtained in this work with 
experimental measurements suggest that the model (initially obtained in a mean-field approach) 
needs to be modified. 
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I. INTRODUCTION 



We study the magnetic properties and critical behavior of a quenched random-exchange 
Ising model, through extensive Monte Carlo simulations. Our starting point was the model 
proposed in Ref. Ill, consisting of iron and ruthenium atoms randomly distributed in a body- 
centered cubic (bcc) lattice, with probabilities 1—x and x, respectively. In this model, atoms 
are treated as Ising spins, Fe-Fe bonds are ferromagnetic (FM), with exchange integral J, 
while Fe-Ru and Ru-Ru interactions are antiferromagnetic (AF) with exchange integrals 
— XJ and — £J, respectively, where A and £ depend on the ruthenium concentration, x, as 
follows: A = £ = £o — £,\X- The parameters £o and £1 were determined from the experimental 
values for the critical temperatures of the system for some concentrations of ruthenium, 



reported in Ref. 2], by fitting these data to the mean-field solution of the model. 



The main goal o 
troduced in Ref. [l 



this work is to establish the universality class of the spin model in- 
through extensive numerical calculations of some of its thermody- 
namic quantities. We used the metropolis algorithm to generate the data and employed 
re-weighting techniques and finite-size scaling (FSS) analysis. We also compare the behav- 
ior of thermodynamic quantities obtained by Monte Carlo simulations with experimental 
and mean-field results. We compare the values of T c (x) with the experimental measure- 
ments and theoretical estimates [l] to determine if the model for the Fei-zRur alloys 
remains adequate in a non-mean-field approach. In this later work, the parameters of the 
proposed Hamiltonian were obtained through a fit of the theoretical critical temperatures, 
obtained in a mean-field-like approximation, with the measured values reported in Ref. Q 
The mean-field result is obtained using the Bogoliubov inequality, in a similar fashion as the 
one employed in the study of order-disorder transitions through a cluster variational method 

a. 

An important aspect of this work is the FSS analysis to determine the critical exponents. 
The study of this model is mainly an investigation of the effects of quenched disorder on 
the critical behavior of a system of Ising spins, like so many theoretical and experimental 
examples presented throughout the last three decades or more. An interesting fact of these 
disordered systems is the theoretical prediction that the introduction of disorder must change 
the universality class of the system if a > 0, where a is the specific heat critical exponent 
for the pure system. This result is known as the Harris criterion [4j. The a > condition 



is satisfied for the case of the pure Ising model in three dimensions, so we can expect three- 
dimensional disordered Ising models to belong to another universality class, when compared 
to their uniform counterparts. However, the Harris criterion says nothing about the new 
universality class. Renormalization-group arguments infer that the universality class of the 
disordered three-dimensional Ising model does not depend on the concentration x. Moreover, 
experimental results show that, for small concentrations x of AF impurities (as in Fei-^Mn^, 
Fei_ a; Zn a ;F2 and Fei_ x Rua; alloys) or nonmagnetic impurities (as in Fei-^Al^ alloys), these 
systems have a continuous transition between a ferromagnetic and a paramagnetic phase 
at a temperature T c (x) < T c (x = 0%), with a critical behavior clearly different from the 
case without disorder (x = 0%), but that seems to be independent of the concentration 
x, within the experimental precision [5j. By small concentrations we mean x < x c for 
nonmagnetic impurities, where x c is the critical concentration above which there is only a 
single paramagnetic phase 6). In the case of AF impurities, we mean x < x g , where x g is the 
concentration for which the transition line between paramagnetic and ferromagnetic phases 
meets the transition line between paramagnetic and spin-glass phases, where frustration is 
relevant Q. 

It is curious that, contrary to experimental results, Monte Carlo simulations presented 
a wide range of different values for the critical exponents of these disordered systems, all 
seeming to indicate some sort of dependence between universality class and impurity con- 
centration. Only recently (late 1990's onward), with the growth of processor capacity and 
available resources for simulations, is theoretical research heading toward a solution to this 
apparent inconsistency, and studies have shown that critical exponents are independent of 
impurity concentration along the transition line between paramagnetic and ferromagnetic 
phases in disordered Ising models fsLJlO | . These works have shown the importance of 
a careful analysis of the scaling behavior of thermodynamic functions, taking into account 
finite-size correction terms due to irrelevant fields in the Hamiltonian, for a correct assess- 
ment of the critical exponents of these systems. 

At first glance, one could suspect that our model does not belong to any of the previously 
discussed universality classes. However, we could think of it as a member of a more general 
family of models with bonds +J and —rJ randomly distributed with probabilities 1—p and 
p respectively. We obtain the Ising model with randomly diluted bonds (RBIM) as r = 0, 
the Ising model with ±J random bonds for r = 1 and our model for some combinations of 



r and p. Note also that the Fei-^Ru^ system we study in this work has a crossover to the 
Ising Model with randomly diluted sites (RSIM) when x = 10%. Therefore, we make the 
hypothesis that our model belongs to the RSIM-RBIM universality class, as is the case also 
for the paramagnetic-ferromagnetic transition line of the ±J Ising model jf], loj]. In this 
sense, our Fei_ x Ru x system would be somewhere in between the r = and r = 1 limits. 
To assess this possibility, we have performed extensive Monte Carlo simulations, with the 
important inclusion of correction-to-scaling terms, as discussed in the previous paragraph, 
in order to estimate the critical exponents and critical temperatures of the model introduced 
in Ref. Ill This model is explained in Sec. [TTJ the simulation and data analysis methods 
are discussed in Sees. IIHI and HVl respectively and our results are presented and discussed in 
Sec. El 



II. MODEL 

Ref. [2I shows that Fe^^Ru^ alloys are found in the bec structure for x < 30% whereas, for 
x ~ 30%, the system undergoes a crystallographic transition to an hep structure. While in 
the bec structure, the lattice parameter increases steadily with the Ruthenium concentration 
x and the system has a ferromagnetic-paramagnetic phase transition at T c (x), as seen in 
Tab. IB In Ref. Q the authors present a detailed discussion of experimental results and first- 
principle electronic-structure calculations that provide evidence that ferromagnetic Fe-Fe 
bonds and antiferromagnetic Fe-Ru and Ru-Ru bonds should effectively model the Fei-^Ru^ 
system. 

TABLE I. Experimental values of critical temperatures for Fei-^Rur alloys, taken from Ref. 0. 

x T c (K) 
0% 1043 
2% 968(2) 
4% 928(2) 
6% 908(2) 
10% 838(2) 



The model proposed by Paduani and Branco (2008) consists of Fe and Ru atoms randomly 
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distributed on a bcc lattice with probabilities 1 — x and x respectively. Each atom has a 
magnetic degree of freedom that is assumed to behave as an Ising-like spin, so the system 
is described by a spin- 1/2 Ising Hamiltonian 

H = -J2 J ijWv (!) 

where the sum goes over all nearest-neighbor pairs, Oi — ±1 for all sites i, and the exchange 
integral Jij takes the values J for Fe-Fe bonds, — A J for Fe-Ru bonds, and — £J for Ru-Ru 
bonds. 

A mean-field solution [1] of Hamiltonian (JT]), using the Bogoliubov inequality 11, 12], 
yields 

[l-xf 2x(l -x) x 2 14 



1 + exp(-2J/k B T c ) 1 + exp(2XJ/k B T c ) 1 + exp(2£J/k B T c ) J 7' 

In Ref. Q, the authors assume that A = £ = £o — £1^ an d present a least-square fit of Eq. 
(E]) to the experimental values of T c in Tab. [B The values reported for the parameters are 
£0 = 0.54(2) and £1 = 5.4(4), so that the final form of the Hamiltonian reads: 

U = -J^djaidj (3) 

where 

CFeFe = 1 

Cy = <! CFeRu = -(0.54 - 5.4x) (4) 
CruRu = -(0.54 -5.4x), 
and the spin variables take the values ±1 and the sum goes over all nearest-neighbor 
pairs. 



III. MONTE CARLO SIMULATIONS 



within 

to 



J I wit 

3, Q 



We studied the three-dimensional Ising system described by the Hamiltonian (|3 
a Monte Carlo (MC) approach. We have employed the metropolis algorithm 
simulate bcc lattices with 2L 3 sites, periodic boundary conditions, and several system sizes 
L, ranging from 5 to 50. Each site on the lattice is randomly chosen to be an Fe or Ru 
atom with probabilities 1 — x and x, respectively. All random numbers were generated using 
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16] . 



ausworthe (shift-register) generator [14| with "magic numbers" p = 1279 and q = 1063 



A Monte Carlo step per spin (MCS) corresponds to N = 2L 3 attempts to flip a single 
spin, chosen at random, or a full sweep on the lattice, when we perform sequential updates. 
As the sequential algorithm proved to be far more efficient in generating independent states 
than the more traditional random choice, it was the one used for most of our simulations. 
The comparison between these two procedures is presented in Subsec. IV Bl 

Our largest simulations ran up to 4.5 x 10 7 MCSs and we made sure to generate at least 
n = 8000 uncorrelated states, with n given by 

n = nMC ^~H (5) 

where umcs is the number of MCSs, t eq is the equilibration time and r is the largest correla- 
tion time estimated by the integral method [3]. As it is the case for the metropolis dynamics, 
the largest correlation time is, in general, obtained from the magnetization autocorrelation 
function. Therefore, in our case, r refers to the magnetization correlation time. 

For practical reasons concerning data storage, we have only saved simulation data every 
10 MCS. However, we have rescaled our dynamic calculations and present all correlation 
times in units of one MCS. This rescaling introduces a systematic error on the integral 
correlation time when typical values of r are close to 10 MCS or smaller, as shown in Sec. 
IV Bl This error is due to the numerical integration using the trapezoid rule, as the integral 
of an exponential decay approaches the exact value from above with numerical error of the 
order ~ 0(n~ 2 ) [l7|. where n is the number of bins. We checked for systematic deviation by 
comparing pairs of r estimates obtained from the same simulation considering a time unit of 
1 MCS and 10 MCS. We found no systematic error for L > 15 and, since r is overestimated, 
it does not affect the precision of our equilibrium analysis. 

When performing Monte Carlo simulations it is customary to work with the dimensionless 
coupling constant K = J/{k B T) and dimensionless temperature T = 1/K. We also define 
other dimensionless thermodynamic quantities such as the extensive energy 

magnetization per spin 



m =^ZX ( 7 ) 



1=1 



specific heat 



magnetic susceptibility per spin 



the quadratic cumulants 



and 



c = ^ [(E 2 ) - (E) 2 } , (8) 



X = f [(m 2 ) - (\m\) 2 } , (9) 



(10) 



[<m 2 >] 

and the more traditional Binder's cumulant 

U = 1 - ^4. (12) 

On all equations above, (• • • ) denotes thermal averages whereas [• • • ] denotes averages over 
disorder configurations. 

For each pair of x and L values, we ran simulations at different temperatures near the 
critical point (up to 20 temperatures for smaller lattices and 5 for L > 30) and used the 



multiple-histogram method [7, 



18 



19] to compute quantities of interest over an almost contin- 



uous range of temperatures. The thermal error associated with those quantities is estimated 
by dividing the data from each simulation in blocks and repeating the multiple-histogram 
process for each block. The errors are the standard deviation of the values obtained for a 
given quantity on the different blocks. For each set of parameters (L, x, T) we average over 
Ns samples of atomic disorder. We chose 10 < N$ < 20, such that the error due to disorder 
was of the same magnitude or smaller than the thermal error obtained for each disorder 
configuration. Finally, we sum both thermal and disorder errors for an estimate of the total 
error. 



IV. DATA ANALYSIS 

In MC simulations we necessarily deal with finite systems. However, our interest lies 
in critical phenomena, which happen in the thermodynamic limit. The critical behavior of 
such systems may be extracted from finite systems by examining the size dependence of the 
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singular part of the free energy density [12[ • In this finite-size scaling approach we write the 
free energy density for a system of linear size L near the critical point as 

f s[ng (t,h,L) ~ L-tfittf^hL^^UiL-^}) (13) 

where t is the reduced temperature and is given by (T — T c )/T c , H is the external magnetic 
field, and h = H/(ksT). We assume that t and h are the only relevant fields while Ui are 
irrelevant perturbations, such that Ui > 0, which ensures that in the thermodynamic limit, 
as UiL~ u)i — > 0, our free energy is a function only of relevant fields. 

Taking appropriate derivatives of the free energy, it is possible to show that some ther- 
modynamic quantities Q (such as magnetization, specific heat, and magnetic susceptibility) 
may be written in the following scaling form: 



Q = Q L e ll+J2Qi L ~ u>i + --^ 



(14) 

where 9 is related to the traditional exponents. As (fT4"j) will be used to fit numerical data to 
estimate critical exponents, we have to truncate the sum at some point. As each additional 
exponent Ui taken into account will add two free parameters to an eventual fit and that 
reduces drastically the number of degrees of freedom, we considered only the first exponent 
oj\ = to. We then have 

m = m L- p/u { 1 + miL-"} , (15) 

X = Xo^ /v {l + XiL- u } > (16) 
c = c L a / u {l + Cl L-"}. (17) 

A similar scaling behavior near the critical point holds for the derivatives 

^ = G L 1 ' V {1 + G 1 L-*'}, (18) 

where G stands for the Binder's cumulant U or quantities linked to the magnetization, such 



as In (|ra|), In (m 2 ), and In (|m| n ) 15]. Equation (I18j) is particularly useful to determine the 



exponent v. Finally, for the critical temperature we have 

T C (L) =T C + AqL" 1 ^ {1 + A X L~"} , (19) 

where T C (L) is the pseudo-critical temperature for a given system size L, and T c is the true 
critical temperature. 
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We perform least square fits using expressions (|T5l) - (|T9l) with all quantities computed 
at several temperatures close to T c . These estimates of T C (L) are obtained locating the 
temperatures T m where we find the maxima of several quantities which diverge as L — > oo (c, 
X, jig) as well as the temperatures Tf at which the cumulants U±, U22, and Ud = U4—U22 
assume particular fixed values. The fixed values used were MC estimates of the universal 
values these cumulants assume at the critical point on randomly site-diluted (RSIM) and 
jond-diluted (RBIM) Ising models: U A = 1.648(3), U 22 = 0.148(1) and U d = 1.500(1) 
6). The calculation of the relevant quantities for different temperatures is done using the 
multiple-histogram method. Examples of this procedure are shown in Fig. [TJfor the specific 
heat and for the derivative of the Binder cumulant, both used to locate different estimates 
of T m . The comparison between the results from both T m and Tf FSS methods provides an 
additional way to check if our model belongs to the RSIM-RBIM universality class. 



Equations (|T5l) - (|T9|) all have four free parameters to be adjusted in the fitting process and 
no stable fits were possible for our data, meaning our statistics should be increased if we 
desire to obtain the exponents a, (3, 7, z/, and uj independently. Since we are more interested 
in obtaining a, /3, 7, and v than in finding precise correction-to-scaling exponents u, we 
employ a procedure similar to the one presented in Ref. 



15 



in which we set a fixed value for 



exponent u and perform a fit with three free parameters, instead of four. Than we change 
the fixed value of u and keep performing fits to obtain the values of the other exponents. 
Once several fits are made, we locate the interval of values of u such that we minimize the 
values of the variance of residuals, x 2 /DOF, where DOF is the number of degrees of freedom 
of the fit. We repeat this procedure using system sizes L min < L < 50, with L min = 5, 10, 
12, 15, 18, 20, and 25 to obtain the L min that globally minimizes x 2 /DOF. Once the best 
L min and corresponding series of values for uj are located, we average the values obtained 
for the exponents, making sure to include only fits with x 2 /DOF < 1.0 in this statistical 
analysis. 



V. RESULTS AND DISCUSSION 



A. Simulation strategy and preliminary results 

For the pure case (x = 0%), we ran simulations at T = T C HTS = 6.35435, which corre 
sponds to the high-temperature series estimate A^ TS = 0.1573725(6) of the critical couplin. 



Img 



for the pure Ising model on a bcc lattice [20|. We used the single-histogram method [7|, |18| 
to locate the temperatures where the peaks of the thermodynamic quantities of interest oc- 
curred. From the peak locations we chose temperatures to perform new simulations. Finally, 
we employ the multiple-histogram re-weighting using the data from at least five simulations 
at different temperatures such that all peaks were found within the interval between the 
minimum and maximum of those temperatures. 

A similar procedure was employed for the disordered case (x = 4%, 6% and 8%). However, 
since no previous estimates for T c were available, we performed test simulations over a wide 
range of temperatures, in order to make a rough estimate of the location of the critical 
region. Figures |2] and |3] exemplify this initial attempt to determine T C (L) for x = 6% and 
L = 30. Once determined a suitable temperature interval for each concentration x and size 
L, we divided it in five (L > 30) to eleven (smaller lattices) temperatures, to simulate and 
employ the multiple-histogram method as discussed above. 

One interesting aspect of Fig. El besides the additional T c estimate it provides, is the 
behavior of the magnetization as a function of temperature. We observe a slight decrease in 
the total magnetization at low temperatures, when we approach T = 0. This decrease is also 
present in the mean-field approximation to this model (see Fig. 2 in Ref. Although not 
as pronounced as the effect presented in Ref. Q, it also occurs in our Monte Carlo approach; 
thus, it is not a mean-field-only effect. For comparison between Fig. [3] and Fig. 2 in Ref. 
[l, it is important to stress that the lowest temperature we simulated is T = 0.5, which 
corresponds to T w 80 K (this value is obtained using J = 14.16 meV, as discussed in Sec. 

EH- 

It is possible to present a qualitative description of this effect. The probability that a 
Ru atom is completely surrounded by Fe atoms as nearest-neighbors is (1 — x) 8 , which is 
quite high for low Ru concentrations (« 51% for x = 8%, ~ 61% for x = 6%, ~ 72% 
for x = 4%, w 85% for x = 2% etc). Next, let us consider a scenario in which the great 
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majority of Ru atoms have no Ru first neighbors. In a situation like this, there is close 
to no frustration and for T = we expect almost all Fe spins to point in one direction 
while almost all Ru spins point the other way. This results in a total magnetization per 
spin smaller than 1 at T = 0. In fact, if we could guarantee that no Ru atom has a Ru 
neighbor, we would have = 1 — 2x. As the temperature rises, thermal fluctuations 

cause spins to reverse. Now let us compare two situations: an Fe atom surrounded by Fe 
atoms, such that all spins are aligned, and a Ru atom surrounded by Fe atoms, such that the 
Ru spin is in the opposite direction to the Fe spins. The probability to reverse the Fe spin is 
proportional to exp [SJ/ksT] while the probability to reverse the Ru spin is proportional 
to exp [— 8A J/k B T}. As our typical values of A are between 0.0 and 0.54, we expect the 
Ru spins to start reversing more easily than the Fe spins; thus, the value of grows 
with increasing temperature. However, at higher temperatures the magnetization should 
decrease, since thermal fluctuations are then strong enough to flip Fe spins surrounded by 
Fe atoms. 



B. Metropoplis dynamics: choosing spins randomly or sequentially 



In the metropolis algorithm, as originally proposed 13| . the choice of the atom to be 



tested for a possible change is usually random. There are, however, other ways to generate 



the Markov chain. It can be done sequentially or, as for the multispin coding version |2l), 
many spins in a sublattice are tested simultaneously. The above methods are all ergodic 
and satisfy detailed balance pj so we ought to choose the most efficient one. 



Multispin coding leads to a drastic reduction in computational time for Ising systems [21 1 
but it is only applicable when it is possible to express the energy difference between any two 
given states as an integer, which is impracticable for our model, since our exchange integrals 
assume noninteger values. We propose an alternative to this method which consists of 
dividing the lattice in two sublattices, the same way as multispin coding; however, instead 
of testing all spins at once we run over the first sublattice testing all possible spin flips 
sequentially and then go to the second sublattice and repeat the procedure. From now 
on we will refer to this update scheme as sequential update metropolis, as opposed to the 
traditional random update metropolis. 

To verify if this sequential method is efficient, compared to the standard random-update 
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metropolis, we obtain a rough estimate of the dynamical exponent z for both methods. This 
is done by fitting our data, at the critical temperature, to the expression 



t = AL Z , (20) 

where A is a constant. This is how the correlation time r is expected to behave at T c , for 
sufficiently large system sizes L. 

In Fig. HI we present fits for the correlation times. We note that the curves for both 
random and sequential correlation times have almost the same slope, which indicates that 
random and sequential algorithms have approximately the same dynamic exponent z, as 
expected. However, random correlation times are always larger than sequential ones for a 
fixed system size. In fact we also plot the ratio 

R T = ^£ (21) 

^seq 

as a function of L and find an almost-constant R T ~ 4, as shown in the insets of Figs. H] (a) 
and H] (b). For L > 15, we have fit the data to R T = A + BL and obtained A = 4.11(8), 
B = 0.000(3) for the pure case and A = 3.8(2), B = 0.003(6) for the disordered case. 

The estimates we get for the dynamical exponent for the pure case are z = 2.02(1) for 
random updates and z = 2.01(1) for sequential updates. For the disordered case we have 
only estimated z for x = 6%, which gave us the figures z = 2.06(3) and z = 2.02(2) for 
the random and sequential update dynamics, respectively. All values are in agreement with 
the MC result, z = 2.04(3), for the three-dimensional (3D) pure Ising model, presented in 



Ref. 
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So, within error bars, both pure and disordered systems have the same dynamic 
exponent for random or sequential updates. However, we find that the correlation times 
for a given system size L are always greater for the random update scheme. Therefore, the 
sequential update version is more efficient than the random update version. 



C. Critical exponents 

Figure shows the comparison between the behavior of some quantities dG/dK at differ- 
ent estimates of the critical point, obtained by the methods described in Sec. [TV] Note how 
the quantities assume almost the same value for the same size at the four different T c esti- 
mates. Therefore, the four different fits for each quantity are almost indistinguishable. As 
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a result, the numerical values we obtain for each critical exponent by independent methods 
are the same, within error bars. Thus, we combine both T m and Tf FSS methods to obtain 
our final estimates of a, 7, (3, and v for the disordered systems. 



TABLE II. Estimates of 1/u using Eq. (|18p for some quantities G, with -L m i n = 15. The region 
0.30 < lj < 0.36 minimized the values of x 2 /DOF. In this interval we found no difference between 
X 2 /DOF up to the second decimal figure for each quantity considered. The values obtained were 
X 2 /DOF = 0.22, 0.23, and 0.081 for G = In (\m\), G = In (m 2 ), and G = U, respectively. 



1/u 




Ct = m \|m|y 


Lj = In {m ) 


G = U 


0.300 


1.465(31) 


1.463(29) 


1.444(47) 


0.305 


1.466(31) 


1.464(29) 


1.445(46) 


0.310 


1.467(31) 


1.465(29) 


1.446(46) 


0.315 


1.468(31) 


1.466(29) 


1.448(46) 


0.320 


1.469(31) 


1.467(28) 


1.448(46) 


0.325 


1.470(30) 


1.468(28) 


1.449(45) 


0.330 


1.471(30) 


1.469(28) 


1.450(45) 


0.335 


1.472(30) 


1.470(28) 


1.451(45) 


0.340 


1.473(30) 


1.471(28) 


1.452(45) 


0.345 


1.474(30) 


1.472(28) 


1.453(44) 


0.350 


1.474(30) 


1.473(27) 


1.454(44) 


0.355 


1.475(29) 


1.474(27) 


1.455(44) 


0.360 


1.476(29) 


1.474(27) 


1.456(44) 



To estimate v we fit the data to Eq. (jl~8{) . where we used G = U, G = ln(|m|), and 
G = ln(m 2 ). For x = 0% we computed the derivatives at the temperature T m where we 
found the maximum of each quantity. Following the procedure described in Sec. IIV} we 
found good fits with L min = 12 for the logarithmic derivatives and L min = 5 for the Binder's 
cumulant. The fits in Fig. 0(a) were done with u ~ 1.0, consistent with the result reported 



in Ref. Il5|. The value we obtain in this work, v = 0.6269(20), is in excellent aer< 

nnnnr 

others in the literature for the three-dimensional Isine; model [5|, ll5|, l2fj|, l23H25[ . 



lent agreement with 
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For the disordered cases, fits using Eq. (|T8|) were done with the same quantities G = U, 
G = In and G = In (m 2 ). Tab. HT1 shows our results for l/v for x — 6 %, obtained from 

data at temperatures T m . One of these fits is shown in Fig. E] (b). Combining the various 
estimates in Tab. HI we find l/v = 1.4632(90) or v = 0.6834(42). Using both T m and T f 
methods, we arrive at our final estimate, v = 0.6826(46), presented in Tab. II III along with 
our final estimates of v for the other Ru concentrations. 

Following the same procedures, we used Eq. f|T7|) to find a. For x = 0%, the results 
are presented in Tab. II I II and compared to other results in the literature. Figure [7J (a) 
shows one of the fits for the specific heat of the pure system: we obtain the estimate 
ajv = 0.1743(70). For the disordered systems, however, no stable fits were possible with 
only one correction-to-scaling exponent and we lack statistical resolution to perform fits with 
higher-order correction-to-scaling terms. Fig. [7] (b) shows one of the plots for the specific 
heat of the disordered system (x = 6%) in which the dashed line is only a guide to the eye. 

The estimates for a for the disordered case, presented in Tab. IHIt were not determined 
independently, but were obtained using the Josephson equality a = 2 — ud, where d is the 
dimension (d = 3 in our case). The and 7 values were all obtained independently by 
fitting Eqs. ( fl~5l) and (fi~6l) . respectively, and our final estimates are presented in Tab. IIHI 
We see that, within error bars, critical exponents do not depend on the concentration x, 
as expected. Note that the usual scaling relations for the exponents are satisfied, for both 
pure and disordered cases. This result is an independent check of our calculation. Also, our 
results for the critical exponents agree, within error bars, with values reported previously 
elsewhere, obtained both from theoretical methods or from experimental studies. This gives 
support to our hypothesis concerning the universality class of the model treated in this work. 



D. Critical temperatures 

For x = 0%, we fit Eq. f|T9|) for the temperatures where we located peaks of several 
thermodynamic quantities. Figure M shows the values of T m used in our FSS analysis. Our 
best fits were obtained with L min = 25 for the magnetic susceptibility and the derivatives of 
(\m\) and U and I/ m ; n = 20 for the remaining quantities. Combining the estimates obtained 
in those fits, presented in Tab. HVt we obtain: T c = 6.3544(6), or K c = 0.15737(1), which 
is in excellent agreement with -^^ TS = 0.1573725(6) [2^] and with another Monte Carlo 
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TABLE III. Our estimates of the exponents a, /3, 7, and v compared to other works in the literature. 
The lines with the labels x = 0%, 4%, 6%, and 8% correspond to our final estimates. The table 
also contains other Monte Carlo [15I.I23I]. 
theory 



24 



high-temperature series expansion 



20| | , and quantum field 



the RSTM JrBTM 



251 ] results for the pure 3D Ising universality class, as well as Monte Carlo results for 



3,3 



n 



5 



and ±J Ising Model [10]. Reference [5j reports experimental results for 
both pure and disordered (with low disorder) 3D Ising universality classes. 









Pure 










a 


P 


7 


V 


a + 2/3 + 7 


2 - a - 3u 


x = 0% 


0.1093(47) 


0.3316(86) 


1.2448(70) 


0.6269(20) 


2.017(29) 


0.010(11) 


[151 


0.1190(60) 


0.3258(44) 


1.2390(71) 


0.627(2) 


2.020(21) 


EE (£ 


[231 


0.1099(7) 


0.32648(18) 


1.2371(4) 


0.63002(23) 


= 2^ 


= 


[201 


0.1094(45) 


0.3266(10) 


1.2375(6) 


0.6302(4) 


EE 2 


EE 


[24, 251 


0.1100(45) 


0.3270(15) 


1.2390(25) 


0.6300(15) 


2.003(10) 


= 


[5] 


0.110(5) 


0.325(5) 


1.25(2) 


0.64(1) 


2.01(4) 


-0.03(4) 








Disordered 










a 


P 


7 


V 


a + 2/3 + 7 


2 - a - 3u 


x = 4% 


-0.062(25) 


0.3584(96) 


1.334(20) 


0.6873(84) 


1.989(64) 


= 


x = 6% 


-0.049(15) 


0.359(16) 


1.330(18) 


0.6826(46) 


1.999(65) 


EE 


x = 8% 


-0.057(20) 


0.3588(81) 


1.337(21) 


0.6856(66) 


2.003(46) 


EE 


[8] 


-0.051(16) 


0.3546(28) 


1.342(10) 


0.6837(53) 


1.994(32) 


0.006(32) 


[9] 


-0.049(9) 


0.3535(17) 


1.342(6) 


0.683(3) 


EE 2 


EE 


[6] 


-0.049(6) 


0.354(1) 


1.341(4) 


0.683(2) 


EE 2 


EE 


[10] 


-0.046(6) 


0.329(2) 


1.339(7) 


0.682(2) 


1.951(17) 


EE 


L5] 


-0.10(2) 


0.350(9) 


1.31(3) 


0.69(1) 


1.91(7) 


0.03(5) 



a The label = corresponds to cases where a was calculated using the Josephson equality. 

b The label = 2 corresponds to cases where either f3 or 7 was calculated using the Rushbrooke equality. 



evaluation [26[, K c = 0.157371(1). We can use our T c to estimate the exchange integral 
of the Fe-Fe bond, JFeFe = J, just by combining the experimental T c in Tab. |T] with the 
definition T = ksT/J. We obtain J = 14.16 meV, which is close to 12.9 meV, reported in 
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Refill and in the interval 10 to 50 meV, as expected for Fe, Co, and Ni [271 ] . 



TABLE IV. Estimates of T c (x = 0%) obtained with FSS analysis of T C (L) for some thermodynamic 
quantities. 



Quantity L n 



c 20 6.35514(34) 

X 25 6.35467(11) 

H 25 6.35447(37) 

25 6.35465(34) 

^ln(|m|) 20 6.35360(57) 

^p- 20 6.35437(26) 

^ln(m 2 ) 20 6.35359(58) 



TABLE V. Our Monte Carlo estimates of critical temperatures compared to experimental and 
mean-field results. 

This work Experimental [2] Mean-field [1] 



x f c T C (K) T C {K) T C {K) 

0% 6.3544(6) - 1043 1043 

2% - - 968(2) 983 

4% 6.1089(23) 1002.7(4) 928(2) 933 

6% 5.9615(17) 978.5(3) 908(2) 893 

8% 5.8064(22) 953.0(4) - 863 

10% - - 838(2) 842 



Following the same procedure for the disordered systems and also including the temper- 
atures Tf (fixed values of the cumulants) in the analysis, we find the critical temperatures 
for x = 4%, 6%, and 8%. The value J = 14.16 meV is used to calculate all critical tempera- 
tures in Kelvin, in order to compare with experimental results, as presented in Tab. |V] The 
mean- field critical temperatures are easily computed with Eq. (T5]), which is taken from Ref. 
Q [see Eq. (14) in Ref. Q. 

Our Monte Carlo estimates of T c for the disordered cases do not agree with the mean-field 
prediction; neither do they agree with the experimental values. We ascribe this discrepancy 
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to the method by which the parameters of the model were obtained. Although those param- 
eters fit the experimental data well in the mean-field approach, the model is not suitable for 
these Fei_ x Rua; alloys in the context of Monte Carlo simulations. 



VI. CONCLUSION 

In this study we used Monte Carlo simulations to investigate the magnetic properties and 
critical behavior of a model proposed by Paduani and Branco (2008) for Fei_ x Ru x alloys 
through a mean-field approach. Our simulations were restricted to ruthenium concentrations 
x = 0%, 4%, 6%, and 8%. We employed re-weighting single and multiple histogram meth- 
ods and finite-size scaling analysis, considering up to the first-order correction-to-scaling 
exponent in order to obtain the critical temperature and critical exponents of the model. 

In the pure case, x = 0%, the values obtained for the critical exponents are in excellent 
agreement with the ones in the literature. The critical temperature we found for the pure 
system also agrees very well with the high-temperature series expansion estimate by Butera 
and Comi (2000) and with another Monte Carlo result [26]. In the disordered cases, we 
show that the critical exponents are consistent with the universality class of the transition 
line between paramagnetic and ferromagnetic phases of three-dimensional Ising models with 
random site or random bond dilution, as well as with the three-dimensional Ising model 
with randomly distributed +J and — J exchange integrals. 

For x = 4%, 6%, and 8%, our estimates of T c do not agree with the mean-field prediction 
nor with experimental results. This is expected, since the model parameters were determined 
through a fitting procedure of experimental data in a mean- field approach However, the 
model proposed in Ref. Ill is in the universality class of three-dimensional disordered Ising 
models. Furthermore, our values for the critical exponents do not depend on x, as expected. 
This result is obtained only if we take into account a correction-to-scaling term in the finite- 
size-scaling analysis. 

To propose a model that brings simulations and experimental results closer together, 
we must seek other ways to determine the dependence of the exchange integrals of Fe-Ru 
and Ru-Ru bonds with the ruthenium concentration. One possibility would be a mean-field 
renormalization group approach, which is presently being carried out. 
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FIG. 1. Thermodynamic quantities versus dimensionless temperature T for x = 0% and lattice 
sizes ranging from 20 to 50. In (a) we plot the specific heat and in (b) the slope dU/dK of the 
Binder's cumulant. The lines were obtained with the multiple histogram method, the symbols 
represent the location of the peaks for each lattice size and the dashed line connecting the peaks is 
a guide to the eye. The tick in T = 6.3544 corresponds to our estimate of the critical temperature. 
The data in this figure and in all other figures in this work were plotted with their respective error 
bars; however, some error bars are smaller then the symbols. 
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FIG. 2. Specific heat versus T for x = 6% and L = 30. The points represent the simulation data 
and the dashed line is a guide to the eye. 
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FIG. 3. (Color online) Total magnetization per spin and individual contributions of Fe and Ru 
atoms for the magnetization versus T, obtained for x = 6% and L = 30. The points represent the 
simulation data and the lines are only a guide to the eye. 
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FIG. 4. (Color online) Log-log plots of the magnetization correlation times versus L for (a) x = 0% 
and (b) x = 6%. The dashed lines are fits using Eq. (|20p . Each inset is a plot of the ratio between 
random and sequential correlation times versus L for the respective ruthenium concentration where 
the solid line is a fit using R T = A + BL. All fits were done for L > 15, to avoid the systematic 
error introduced by numerical integration, present for correlation times smaller than or equal to 
10 MCS, as discussed in Sec. IIIII 
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FIG. 5. (Color online) Log-log plot of the derivative dG/dK versus L for x = 4%, for (a) G = 
ln(|m|) and (b) G = ln(m 2 ) [see Eq. (|18p and remarks after it]. The data were computed at 
slightly different critical point estimates, given by the temperatures T m (maximum of dG/dK) and 
Tf, in order to compare the estimates. The solid lines are fits to dG/dK = a%lM v . The inset 
corresponds to the same graph made with both data and lines slightly shifted along the L axis to 
make them visible. 
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FIG. 6. (Color online) Log-log plots of the maxima of some thermodynamic quantities G as 
functions of L that were used to determine v for (a) x = 0% and (b) x = 6%. The full lines are 
fits performed with Eq. (fT8|) for L m \ n < L < 50. The dotted lines are extrapolations of the fits for 
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FIG. 7. (Color online) Maxima of the specific heat versus L. For x = 0% (a), the plot is in log-log 
scale and the dashed line is a fit of the data using Eq. (fT7|) . with to = 1.0. For x = 6% (b), the 
plot is on a linear scale and the inset corresponds to the same data on a log-log scale. The dashed 
line in (b) is only a guide to the eye. 
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FIG. 8. Effective critical temperatures T C (L) versus L l / y estimated for several quantities. The 
dotted and dashed lines are fits performed with Eq. [T9l and using our estimate v = 0.6269. 
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